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Abstract 

We present a Markov-chain Monte Carlo algorithm of worm type that correctly simulates the 
0{n) loop model on any (finite and connected) bipartite cubic graph, for any real n > 0, and 
any edge weight, including the fully-packed limit of infinite edge weight. Furthermore, we 
prove rigorously that the algorithm is ergodic and has the correct stationary distribution. We 
emphasize that by using known exact mappings when n = 2, this algorithm can be used to 
simulate a number of zero-temperature Potts antiferromagnets for which the Wang-Swendsen- 
Kotecky cluster algorithm is non-ergodic, including the 3-state model on the kagome lattice 
and the 4-state model on the triangular lattice. We then use this worm algorithm to perform a 
systematic study of the honeycomb-lattice loop model as a function of « < 2, on the critical 
line and in the densely-packed and fully-packed phases. By comparing our numerical results 
with Coulomb gas theory, we identify a set of exact expressions for scaling exponents governing 
some fundamental geometric and dynamic observables. In particular, we show that for all n < 2, 
the scaling of a certain return time in the worm dynamics is governed by the magnetic dimension 
of the loop model, thus providing a concrete dynamical interpretation of this exponent. The case 
n > 2 is also considered, and we confirm the existence of a phase transition in the 3-state Potts 
universality class that was recently observed via numerical transfer matrix calculations. 

Keywords: Monte Carlo; Worm algorithm; Loop model 
PAC5.02.70.Tt,05.10.Ln,64.60.De,64.60.F- 



1. Introduction 

Among the myriad of models studied in the theory of critical phenomena, two fundamental 
examples that continue to play a central role are the ^-state Potts model , and the 0{n) spin 

model |4, 5]. In the original spin representation, the parameter qor n must be a positive integer 
However, the Fortuin-Kasteleyn representation |6] of the ferromagnetic Potts model and the loop 
representation [7J of the 0{n) spin model show how these models can be extended to arbitrary 
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real q,n > 0, by re-expressing them as models of random geometric objects: clusters or loops, 
respectively. In fact, the extension of the Potts model to non-integer q can also be formulated 
directly in the spin language, by re-expressing the Potts spin clusters in terms of domain walls (ST . 
These geometric models play a major role in recent developments of conformal field theory ||9|] 
via their connection with Schramm-Loewner evolution (SLE) |10, 11]. 

Monte Carlo methods are an indispensable tool in statistical mechanics 1 13 131 The Sweeny 
algorithm fl4] and the Swendsen-Wang-Chayes-Machta cluster-algorithm llisl Il6ll provide re- 
markably efficient [17^ JJJ tools for studying the ferromagnetic Potts (random-cluster U9n ) 
model, and are valid for any real q > 0, or q > 1, respectively. For loop models, by con- 
trast, efficient simulation at noninteger n has posed a significant challenge. Instead, numerical 
transfer-matrix techniques have typically been used 1,20. 211. Monte Carlo simulations at n 1 
have been reported in 11221 12311 . however the algorithms used were essentially single-spin-flip 
Metropolis algorithms for Ising spins on the dual lattice. As such, their efficiency is limited, 
and they are manifestly non-erg odicQ at infinite edge weight. In |24], a cluster algorithm was 
presented that is valid for all n > 1, which is dramatically more efficient than the single-spin-flip 
algorithms on the critical line. However, its efficiency deteriorates rapidly as the edge weight 
increases, and it too becomes non-ergodic at infinite edge weight. 

In II25II . a Monte Carlo algorithm of worm type for simulating the honeycomb-lattice fully- 
packed loop model with n - \ was presented, and its validity was rigorously proved. In this 
article, we present a worm algorithm that correctly simulates the 0{n) loop model on any bipartite 
cubic gr aphll for any real n > 0, and any edge weight, including the fully-packed limit of 
infinite edge weight. Furthermore, we prove rigorously that the algorithm is ergodic and has the 
correct stationary distribution. We then use this algorithm to perform a systematic study of the 
honeycomb-lattice loop model as a function of n. By comparing our numerical results for n < 2 
with Coulomb-gas theory |26], we identify the exact scaling exponents of some fundamental 
geometric observables, as well as certain observables related to dual Ising spins. Furthermore, 
we find that for all « < 2, the scaling dimension of a certain very natural return time in the worm 
dynamics coincides precisely with the magnetic dimension of the loop model, which provides 
a concrete dynamical interpretation of this exponent which is meaningful for both integer and 
non-integer n. See section [3] We also study the case n > 2, and confirm the existence of 
a critical transition in the 3-state Potts universality class, which was recently observed using 
transfer matrices |21]. While the honeycomb-lattice model is perhaps the archetypal loop model, 
and is certainly the most well-studied case, there are other examples of bipartite cubic graphs 
which are of interest, including the (4 ■ 8^) Archimedean lattice (dual of the Union Jack lattice), 
and the Hydrogen-peroxide lattice (which is three-dimensional). Systematic studies of the loop 
models on both of these lattices can be performed using the algorithms described in this article; 
the results will be presented elsewhere. 

In addition to the study of loop models, the worm algorithms that we present here can also 
be appHed to the study of a number of antiferromagnetic Potts models. It is well known that 
the honeycomb-lattice fully-packed loop model with n = 1 is equivalent to the zero-temperature 
triangular-lattice antiferromagnetic Ising model. The latter model (which is critical) provides a 
canonical example of geometric frustration, but is notoriously difficult to simulate. In fact, even 



Following the typical usage in the physics literature, we take ergodic as synonymous with irreducible. Recall that a 
Markov chain is iiTeducible if for each pair of states and j there is a positive probability that starting in i we eventually 
visit j, and vice versa. 

^All graphs considered in this article are implicitly assumed to be finite and connected. 
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the most sophisticated tailor-made cluster algorithms II27L 1281] are thought to be non-ergodic. 
However, the worm algorithm constructed in [25J immediately provides a provably ergodic 
Monte Carlo method for this problem. Furthermore, it is known that the honeycomb-lattice 
fully-packed loop model with « = 2 is equivalent to both the zero-temperature kagome-lattice 
3-state Potts antiferromagnet and the zero-temperature triangular-lattice 4-state Potts antiferro- 
magnet ll29r . Both of these models are believed to be critical. While the Wang-Swendsen- 
Kotecky |30] (WSK) cluster algorithm is undoubtedly the current state-of-the-art for simulating 
antiferromagnetic Potts models, it has recently been proved ifsil [32I1 to be non-ergodic for both 
of these cases. By contrast, the worm algorithms described in Section |2] have been proved to be 
ergodic and can be applied in a straightforward way to the study of both of these Potts antifer- 
romagnets. The details of this application will be reported elsewhere (but see also Section|4]for 
further discussion). 

Worm algorithms were first applied to classical lattice models in f33^, and it was demon- 
strated empirically in Ii34i1 that the worm algorithm is an extraordinarily efficient method for 



simulating the three-dimensional Ising model. See 11351 13611 for some recent applications to 0(n) 
models. Worm algorithms provide a natural way to simulate cycle-space models. Given a fi- 
nite graph G = {V,E), the cycle space, C(G), is the set of all A c £ such that every vertex in 
the spanning subgraph {V,A) is even. We call A c E and (V,A) Eulerian whenever A e C(G). 
Fig.[l |{a)| shows a typical configuration on the honeycomb lattice. The essence of the worm idea 
is to enlarge the state space C(G) to include a pair of defects (i.e., vertices of odd degree), and 
then to move these defects via random walk. When the two defects collide, the configuration 
becomes Eulerian once more. A very natural class of cycle-space models is defined for n, x > 
by the probability measure 

0c,«,.(A)'x n-WxI^I, AeC(G), (1) 

where c(A) is the cyclomatic number of the spanning subgraph (V,A). The cyclomatic number 
of a graph is simply the minimum number of edges to remove from it in order to make it cycle- 
free. On graphs G of maximum degree A(G) < 3 therefore, all A e C(G) consist of a collection 
of disjoint cycles, or loops, and c(A) is then simply the number of such loops. Consequently, 
the model ([1]) is typically referred to as the loop model, and bipartite cubic graphs (such as the 
honeycomb lattice) provide a natural setting for its study. 

It is well known ff\ that on any graph G = {V, E) of maximum degree A(G) < 3, the model 
([T]l arises for positive integer n as a loop representation of an n-component spin model, 

Z = Tr ]~[ (1 + « XO-, ■ o-j), (2) 

ijeE 

where cr - (cr', . . . , cr") e K" and Tr denotes normalized integration with respect to any a priori 
measure (■)o on IR" satisfying {cr^crP)^ - da^/n and (cr"')o = {cr"crPcr^)Q - 0. In particular, 
uniform measure on the unit sphere is allowed, as are various /ace-cM^j/c and corner-cubic mea- 
sures liTJ. For n 5t 1, the Boltzmann weight (|2| with spins on a sphere defines a nonstandard 
0(n) spin model, which has positive weights only for |x| < 1/n, but it is, nevertheless, expected 
to belong to the usual 0{n) universality class. 

In the limit x — > -1-00, the support of 0g,«,x reduces to the set of all A e C(G) with maximal |A|. 
On bipartite cubic graphs maxAec(G) 1^1 = so the set of all such/M/ty-pacfet/ configurations is 
simply 

T{G) ;= {A e C(G) : t/„(A) = 2 for all v e V]. (3) 
3 



Figure 1: Typical loop configuration |(a)| and fully-packed loop configuration |(b)| on the honeycomb lattice with periodic 
boundary conditions. Thick lines denote occupied edges, thin lines denote vacant edges. 



(a) Loop configuration 



(b) Fully-packed loop configuration 



Fig- [J |tb)| shows a typical fully-packed configuration on the honeycomb lattice. We note that 
the elements of TiG) are referred to as 2-factor^hy graph theorists |37]. We also remark that 
A e T{G) iff £■ \ A is a dimer covering (perfect matching) of G. Finally, we note that the limiting 
measure is simply 

0G,«(A) := 0G,«,co(A) (X „'(^>, A € T{G). 

A great deal is known about the loop model ([1]) on the honeycomb lattice when n < 2. For 
given n, the model is believed to have three distinct phases: a disordered phase (small x), a 
densely-packed (DP) phase (large finite x), and a fully-packed (FP) phase (infinite x). Further- 
more, the model is exactly solvable ||38 , 39 , 4^, 41 , 42 ] on the curves 
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(4) 



The plus sign in (|4]i corresponds to the critical curve, Xc(n), separating the disordered and 
densely-packed phases ll43ll . The minus sign in (|4]i corresponds to a curve of stable fixed points 
in the densely-packed phase. For all finite x > x^. , the loop model is in the densely-packed phase, 
which is critical ll20ll . The fully -packed model is also critical, however it is known to be in a 
distinct universality class to the densely-packed phase 144 l45l 14611 . For convenience, we shall 
refer to the model ([TJ with x^ < x < oo as the densely-packed loop (DPL) model, and to the 
X = CO model as the fully-packed loop (FPL) model. 

The loop model with loop fugacity n can be related to a Coulomb gas ll26l. l46ll with coupling 
^by 

n - -2 cos{ng/4), (5) 



with 



8^ 



([2,4], 
I [4, 6], 



Xc < X < +00, 



(6) 



'More precisely, if A g T(G) then the spanning subgraph {V,A) is a 2-factor. 
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Recall II47II that the critical (densely-packed) loop model with < « < 2 corresponds to a 
tricritical (critical) Potts model with q - 11^. The normalization for g given in (|5]l and (|6]l, which 
is a factor of 4 times larger than the g presented in [26], is in fact the standard normalization for 
the Coulomb gas corresponding to the q - rP' Potts model, rather than the 0{n) loop model. This 
choice facilitates easy translation between loop and Potts exponents, which will prove convenient 
in Section[3] 

Coulomb gas theory flS] predicts a whole spectrum of exact scaling dimensions character- 
izing the loop model. However, we emphasize that identifying which loop-model observables 
these exponents actually govern is not always obvious. In |24|], Monte Carlo simulations were 
combined with Coulomb gas predictions to identify the exact scaling exponents for a number 
of natural geometric observables, as well as observables related to dual Ising spins. The results 
presented in 1241. however, were restricted to the critical branch, x - Xc, and to n > 1. In this 
work we shall use worm algorithms to extend these observations to all n > 0, and to the DP and 
FP phases. See Section[3] 

The outline of this article is as follows. The necessary theoretical results concerning the 
algorithms appear in Section|2] In Section[3]we then present our numerical results for the n < 2 
honeycomb-lattice loop model in the critical, densely-packed, and fully-packed phases, and also 
for the n > 2 model. Section|4]then concludes with a discussion. 



2. Worm dynamics for loop models 

We begin our discussion of worm dynamics by constructing a worm algorithm to simulate 
([TJ on an arbitrary graph, for any < n, x < co. This essentially generalizes the presentation 
in 1I34 , 25 1 to include a loop fugacity in the stationary distribution. We then demonstrate that it is 
possible to make this algorithm rejection-free, in a certain sense. We then turn our attention to the 
specific case of cubic graphs, and consider the limit x ^ 00. In particular, we rigorously prove 
that the rejection-free worm algorithm on any bipartite cubic graph remains ergodic aX x - +00. 



2.1. Simple worm dynamics 

Fix a finite graph G - (V, E), and for any A Q E dA Q V denote the set of all vertices 
which have odd degree in the spanning subgraph (V, A). Loosely, dA is just the set of sites that 
touch an odd number of the bonds in the bond configuration A. If m, v e V are distinct we write 

^„,,.(G) := {A c £ : 5A = {M,y)}, 

and 

5,,„(G) {A c £ : 5A = 0). 

We emphasize that Sy^viG) - C{G) for every v € V. We take the state space of the worm 
algorithm to be 

S{G) {(A, M, v) : M, V e y and A € >S„,„(G)), 

i.e., all ordered triples (A, u, v) with A Q E and m, v e V, such that A e Su^yiG). Note that if 
(A, u, v) € S{G) then A € C(G) iff m = v. Thus, the bond configurations allowed in the state space 
of the worm algorithm constitute a superset of the Eulerian configurations. Finally, we assign 
probabilities to the configurations in S{G) according to 



7!'G,n,x(A, M, v) oc c/„ ^/^, n^*'*' x^\ (A, u, v) € S(G), 
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(7) 



where t/,, denotes the degree in G of v e V. In the following, when we wish to refer to the degree 
of V G y in the spanning subgraph {V,A) we will write dt,(A). Loosely, di,{A) is simply the number 
of bonds that touch v in the bond configuration A. In this notation we have = dv(E). 

The first step in constructing the standard worm algorithm is to consider the worm proposal 
matrix, P^^\ which is defined for all uu' e E and v e V by 

P^°'[(A, M, v) ^ (Aamm', m', v)] = P*">[(A, V, u) {Aauu', v, u')] = — , (8) 

2£/„ 

all other entries being zero. Here A denotes symmetric difference, i.e. delete the bond uu' from 
A if it is present, or insert it if it is absent. It is easy to see that f is an ergodic transition matrix 
on S{G). According to (|8]l, the moves proposed by the worm algorithm are as follows: Pick 
uniformly at random one of the two defects (say, u) and one of the edges emanating from u (say, 
uu'), then move from the current configuration (A, u, v) to the new configuration (A A uu', u', v). 

Now we simply apply the usual Metropolis-Hastings prescription (see e.g. fl3^) to assign 
acceptance probabilities to the moves proposed by P'"', so that the resulting transition matrix, 
Pc,n,x, is in detailed balance with (|7]i. Explicitly, for all uu' e E and v e V we have 

Pc.n.xK^' ^) ^ (AAMm', u', v)] = Pc,n,x[{A, V, u) — > (AAUu' , V, u')] 

min(l,xn) mm' ^ A and m <-> m' in (V, A) 

1 min(l,jic) mm' ^ A and M m' in (y. A) 

2du min(l,l/nx) mm' e A and m «-> m' in (V, A \ mm') 

min(l, l/x) uu' e A and u ^ u' m (V,A \ mm') 



(9) 



The notation m m' in (O means that vertices u and u' are connected in the stated spanning 
subgraph of G. The transitions (|9j define Pc.n. v uniquely since all other transitions occur with 
zero probability, except the identity transitions (A, u, v) — » (A, u, v), whose transition probabilities 
are fixed by normalization. 

Now let us return to our original goal, which was to sample from CiG). As elaborated below 
in Lemma ITTl achieving this is as simple as running the worm chain and choosing to only 
measure observables when the two defects meet, u - v. This defines an ergodic Markov sub- 
chain on 

{(A, M, v) e S{G) : M = v) s CiG) x V (10) 

with stationary distribution 'nc,n,xiA, v) oc n'^*'*' x''^', and therefore for any loop-model observable 
X : C(G) ^ R we have {Xh^J^ = <X)^,.„ , . 

The resulting Monte Carlo algorithm is summarized in Algorithm[T] Note that the acceptance 
probabilities (which are simply 2 du Pc.n.xliA, u, v) — > (Aauu', u', v)]) will in general depend on 
the topology of the loops in a non-trivial way; we shall return to this point in Section 12.61 The 
abbreviation UAR simply means uniformly at random. 

Remark 2.1. The naive n — > limit of ([T) with x < oo held fixed reduces to a trivial measure 
concentrated on the single state A = 0, and the naive n Q limit of Algorithm[T]leads to a trivial 
dynamics which correctly simulates this trivial model. Non-trivial n — > limits can be taken 
however; for example, conditioning on positive cyclomatic number before taking the n — > limit 
of ([TJ yields a model of self-avoiding polygons. Worm dynamics can be developed to simulate 
both self-avoiding walks and self-avoiding polygons, however a discussion of these issues would 
lead us too far afield here. A discussion of such algorithms will be reported elsewhere. 
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Algorithm 1 (Simple worm algorithm). 
loop 

Current state is (A, u, v) 

UAR, pick one of the two defects (say u) 

UAR, pick a neighbor u' ofu 

Make the transition (A, u, v) — > (Aamm', u' , v) with acceptance probability inferred from Q 
if u' - V then 

Measure observables 
end if 
end loop 



2.2. Markov sub-chains 

Let us pause momentarily, and consider, quite generally, that we have an ergodic Markov 
chain on a finite state space iS which is in detailed balance with a distribution n, and suppose that 
we only observe the process when it is in a state in A" c >S. This new process is a Markov chain 
on X with transition matrix 

oo n 
"=0 .so,.!i,...,.!„e^ 

Lemma 2.1. P is ergodic and in detailed balance with the restriction of n to X 

- „ 

TT, - — , i 6 A. 



Proof. The ergodicity of P on A" follows immediately from the ergodicity of P on S, and one 
can easily verify directly that P is in detailed balance with n by using the fact that P is in detailed 
balance with n. □ 

We now wish to make the following observation. Since we only observe the ^S-chain when it 
visits A c iS, we have quite a bit of freedom to modify the transition probabilities in X - S\X 
without affecting the stationary distribution of the A-chain. In particular, we can forbid identity 
transitions in X. 

Corollary 2.2. // 

{iP)ss' seX 



JDlf— seX,s'^s 

1 - (P)s.s 

seX,s' ^ s 



then P' is in detailed balance with n. 

Proof. It is elementary to verify directly that P' is in detailed balance with 

, s^X, 
''■'"\(l-(P),,)7r, seX. 

Lemma |2J] then immediately implies that P' is in detailed balance with n. □ 
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Therefore, both P and P' are in detailed balance with the same distribution n. Since P' 
forbids identity transitions s s when s e X (which correspond to rejections in the context of 
Metropolis algorithms) P' is clearly more efficient than P at sampling from X. We shall refer 
to P' as a rejection-free chain, although it should be emphasized that rejections are still allowed 
inside X. 

A concrete example of the advantage of using the rejection-free chain is provided by consid- 
ering worm algorithms for fully-packed loop models. Indeed, as we shall see, the x — » oo limit 
of Pc,n,x is absorbing, while the corresponding rejection-free algorithm remains ergodic (at least 
on bipartite cubic graphs). 

2.3. Rejection- free worm dynamics 

Thus far we have glossed over an important issue, namely the ergodicity of Pc.n.x- It is not 
hard to see that f G,n,v is ergodic whenever x < oo. However, Pc.n..^ is manifestly non-ergodic 
when X = H-oo, in fact it is absorbing. Indeed, as x ^ oo the probabilities for transitions that 
remove an edge vanish. Consequently, 

D WA ^ Ml du-du(A) d,.-d,.{A) 

f„,+oo[(A,M,v) ^ (A,M,v)] = 1 — — , 

2 du 2 dv 

and aU states (A, u, v) e S{G) for which both du{A) - du and dy{A) - dy become absorbing as 
X —> oa. 

Suppose now that G is A;-regular. Then (A, u, v) will be absorbing when x = +oo iff ii„(A) - 
c/v(A) - k. By definition, if (A, u, v) G S{G) then when u - v the vertex degree t/„(A) - t/,,(A) is 
even, whereas when u + v both ii„(A) and dv{A) are odd. Thus, if k is odd then (A, u, v) can be 
absorbing only \f u v, and so all states with A € C(G) remain non-absorbing. In particular, on 
a cubic graph we can now see that as x — > oo all states (A, v, v) e C(G) x V remain non-absorbing 
while all states (A, u, v) with u v and c/„(A) = dy{A) - 3 become absorbing; once both defects 
have degree 3, the Pcji,x chain remains in the given state for eternity. 

We cannot, therefore, use Pc,n,x to simulate the fully-packed loop model. As we shall see, 
however, we can use its rejection-free counterpart. Following Corollarv 12.21 we define a new 
transition matrix P'q^^ by explicitly conditioning on making a non-trivial transition whenever 
M 4^ V. Specifically, we define 

P'g n ^) ~* (AAMm', U , v)] = P'^j „ ^.[(A, V, u) (AAUu' , V, u')] 

PG,n,x[{A, u, v) (Aauu', u' , y)] U — V, 



f G,n,.v[(A, M, v) (AAUu', u' , v)] 

I 1 - ^'g,«,x[(A, u, v) -> (A, M, v)] 

P'cn v:[(^' ") ~* ")] ~ Pc,n,x[{A, U, u) (A, M, «)]. 



(11) 

U + V, 



All other transitions occur with zero probability. In particular, no identity transitions are allowed 
from non-Eulerian states. Corollary |22] immediately implies that P'^^^ ^ can be used to simulate 
4>c,n,x for any < «, x < oo. We now proceed to show that in fact P'^^ ^ remains valid even at 

X - -Hoo. 

Remark 2.2. While the explicit conditioning (fTTl i is perhaps the simplest way to ensure ergodic- 
ity is retained as x — > oo, there are variations of this idea that also work. Indeed, the algorithm 
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Figure 2: Possible topologies of the defect cluster Cu,-(A) when di,{A) = 3. 




U 



(a) Tadpole graph. 



(b) Dumbbell graph. 



(c) Theta graph. 



presented in 11251] for the n - 1 case was constructed in a sHghtly different way; the main con- 
sequence is that while configurations with t/„(A) = dy{A) = 1 were allowed in ll25ll . they are 
forbidden as x — » oo in the algorithm we present here, as we shall now see. 

2.4. Worm dynamics for fully-packed loops on bipartite cubic graphs 

We now consider the .x — » oo limit of P'q„ ^- We begin by noting that whenever G has max- 
imum degree A(G) < 3 there are only a small number of possible topologies that the connected 
components of states in S{G) can have (see Fig.|2|. 

Proposition 2.3. Let G be a finite graph with A(G) < 3, let (A, u, v) e S(G), and let C„,,(A) be 
the component containing u and v in {V, A). Then all components other than C„y(A) are isolated 
vertices or cycles, and we have the following classification of the possible topologies of the 
component C„i,(A) 



C„,,(A) = 



isolated vertex 


du{A) -- 


-- dM) = 





path 


du(A) -- 


- d,(A) = 


1 


cycle 


duiA) -- 


-- dy(A) = 


2 


tadpole graph 


du(A) -- 


-3,d,(A) 


= 1 or t/„(A) = \,dM) = 3 


dumbbell or theta graph 


du(A) -- 


-- d,(A) = 


3 



Proposition l2.3l is intuitively obvious and its proof (which we omit) is straightforward. 

For the remainder of this section we shall restrict attention to bipartite cubic graphs. On such 
graphs, the limit as jc — > oo of ^ is now easily seen to be given by 

lim f ^[(A, M, v) — > (Aauu , u' , v)] 



1/6 u — V, uu' i. A, 

1/4 t/„(A) = £/,(A) = 1,mm' ^A, 

1/2 t/„(A) = 1, t/,,(A) = 3, uu' i A, 

1 /6 C„,,(A) is a theta graph, 

n/2(n + 2) C,„.(A) is a dumbbell graph, m ^ m' in (V, A \ uu'), 

1 /2(n + 2) C„,.(A) is a dumbbell graph, m <-> m' in (V, A \ uu'), 
9 



(12) 



and 

lim P' J(A, u, u) (A, u, u)] = 2/3. (13) 

X — >co ' '' 

All other transitions occur with zero probability. 
Now consider the subspace 

KiG) := {(A, u, v) e S{G) : dAA) + for all x and c/„(A) + d,{A) > 4). (14) 

If (A, M, y) e %{G), then either |A| = \V\ or |A| = |y| + 1. Furthermore, since (fT2l i only allows the 
deletion of edges when |A| - \V\ + 1, it is clear that 

lim P'q „ J(A, M, v) ^ (Aamm', m', v)] = 

whenever (A,m, v) € %(Cj) and (Aamm',m', v) ^ 'R(G), so "?? is closed (and therefore recurrent). 
The restriction of PJ, ^ to 'R(G) therefore defines a Markov chain on 'R(G). We emphasize that 
the set of all bond configurations A for which (A, v, v) e 'R(G) corresponds precisely with TiG). 

Remark 2.3. Since G is a bipartite cubic graph, we know that |A| - \V\ for all A € TiG). It is 
therefore natural to consider lim^^M x"'^' Zc ji'^ ^(A, u, v), where Zc^n.x is the appropriate nor- 
malization constant (partition function). It is straightforward to verify that this limiting measure 
is supported on fiiG) and is detailed balance with (fT2l i. 

Remark 2.4. The space "RjG) given by ( fT4b is strictly smaller than the state space of the worm 



dynamics considered in II25I1 . In particular, no states in which the defect cluster is a path are 
allowed in ( fT4b . 

Before proceeding further, it is useful to note that ??(G) has the disjoint partition ??(G) - 
SUTU&UD, where 



& = {(A, u, v) e nG) 
T = {(A, M, v) e :R(G) 
= {(A, M, v) e :R(G) 
D = {(A, M, v) € :R(G) 



C„,,(A) is a cycle), 

C„,,(A) is a tadpole graph), ^^^^ 
C„,,(A) is a theta graph), 
C„,,(A) is a dumbbell graph). 



Now let us denote the restriction of (fT2l) to "R(G) by f ^ ^. From (fTSI) we can see that the only tran- 
sitions from (A, M, m) € £ which occur with non-zero probability are (A, u, m) — > (A U mm', m', m) 
and (A, m, m) — » (A U mm', m, m'), which both occur with the same probability 1/6, and the iden- 
tity transition (A, m, m) — » (A, m, m). We are therefore free to multiply the two transition prob- 
abilities P^^[(A,M, m) — > (Aamm', m', m)] = 1/6 from states (A,m, m) e £ by a constant factor 
< y < 3, provided that we also redefine the probabilities for identity transitions, so that we 
retain correctly-normalized row sums. The only effect of such a modification is to multiply the 
stationary probabilities of all the states (A, m, m) e £ by the same constant 1/y. Therefore, such 
modifications do not affect detailed balance. If we now choose 7 = 3, then the identity transitions 
from £ will occur with zero probability. We thus obtain an entirely rejection-free Monte Carlo 
method. 

Putting all these details together, let us now define the following transition matrix Pc,n on 
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Pc,n[{.A., U, v) (AAUU, U , v)] — Pc,n[i.A-, V, u) (AAUU , V, u')] 

'ill (A,M,v) efiUTand mm' ^ A, 

1/6 (A,M,v)e0, 

n/2(n + 2) (A, m, v) e £) and mm' is a bridge, 

1 /2(n + 2) (A, M, v) € D and mm' is not a bridge. 



(16) 



All other transitions are assigned zero probability; in particular, no identity transitions are al- 
lowed. The transition matrix corresponds to a very simple dynamics: if |A| = \V\ we must 
add one of the vacant edges incident to one of the defects; if |A| = |y| + 1 must delete one of the 
occupied edges incident to one of the defects. 

Remark 2.5. Note that the only appearance of n in (fT6] l occurs when du{A) - d^,(A) - 3. Loosely, 
(fTSI l says that for transitions from degree 3 defects, the relative weight given by Pq „ of traversing 
a bridge is n, while the relative weight of traversing a non-bridge is 1 . 

It is now elementary to show that f c,« is in detailed balance with the following distribution 
on "R(G) 

1 (A, M, v) e fi U r, 

3/« (A,M,v)e0, (17) 

in + 2) In (A, m, v) e D. 



7rG,«(A,M, v) 



Z„ 



We prove in the next section that f c „ is also ergodic. It then follows from ( fTTl l and Lemma ITT] 
that Pc,n defines a valid Markov-chain Monte Carlo algorithm to simulate the fully-packed loop 
model on any bipartite cubic graph. We summarize this algorithm in Algorithm|2l 

Algorithm 2 (FPL worm algorithm). 
loop 

Current state is (A, u, v) 
ifu-v then 

Choose the unique edge uu' i A 

Perform, UAR, either (A, m, m) ^ (A U mm', m', m) or (A, m, m) ^ (A U mm', m, m') 
else if u + V then 

if du(A) - 1 or dy{A) - 1 (say u) then 

Choose, UAR, one of the 2 vacant edges incident to u (say uu') 

Make the transition (A, u, v) — > (A U mm', m', v) 
else if du(A) - dy{A) - 3 then 

Choose, UAR, one of the 2 defects (say u) 

Choose with probability P(m') one of the 3 neighbors of u (say u') 
Make the transition (A, m, v) ^ (A \ mm', m', v) 
end if 
end if 
end loop 



We note that the only place in which the topology for the loop configuration enters into 
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Algorithmic] is via the definition of the function P(m'), which from ( fTSb is given by 



1/3 (A,M,v)e0, 

«/(« + 2) (A, M, v) e £) and mm' is a bridge, (18) 
1 /(n + 2) (A, M, v) e £) and mm' is not a bridge. 



See Section 12.61 for a discussion of some possible implementations of the required topological 
queries. 

Remark 2.6. The « — > limit of Algorithm |2] is well defined and non-trivial. Taking the n — » 
limit of ( fT6b or ( fTST i we see that the only change to Algorithm |2] when n = is that we are 
forbidden to delete bridges. If G is hamiltonian, this suggests that the set of recurrent Eulerian 
states should be the set of all Hamiltonian cycles of G. Indeed, the n — > limit of ( fTTb shows 
that 7rG,n=o uniformly samples Hamiltonian cycles, and it is supported on the subset of 'R(G) in 
which the spanning subgraphs (V, A) are connected. 

While it therefore seems plausible that the n — > limit of Algorithm|2]provides a valid Monte 
Carlo method for uniformly sampling Hamiltonian cycles, our general proof of the ergodicity 
of ( fT6b breaks down at « = 0. Although Algorithm |2] is very well suited to simulating fully- 
packed loops for n > 0, it is not perhaps the most natural nor the most efficient dynamics for 
the special case of n = 0, and so we have not attempted to prove its ergodicity in this limit. A 
more natural worm dynamics for simulating Hamiltonian cycles is presented in [48, applying 
earlier ideas from |50]. We expect the specialized dynamics presented in |48, 49] is more efficient 
ffian the n — » limit of Algorithm |2] and it has the added advantage that it simultaneously 
simulates both Hamiltonian paths and Hamiltonian cycles. 

2.5. Ergodicity of Pc,n 

We begin by making some brief comments on our notation. Consider a Markov chain on a 
state space S, with transition matrix P. We say s e S communicates with s' e S, and write 
s-^s' , if the chain may ever visit state s' with positive probability, having started in state s. We 
say states s and s' intercommunicate, and write s'^s' , if s~^s' and s'-^s. Using this notation, 
a set of states J[ Q S is, ergodic iff s'^s' for all s, s' e Jl. Finally, we say P is ergodic if S is 
ergodic under P. 

Proposition 2.4. IfG is a finite, connected, bipartite cubic graph and « > 0, then P^n is ergodic. 

Proof. Let G be a finite, connected, bipartite cubic graph, let B e TiG) be an arbitrary, but fixed, 
fully-packed subgraph, and let z e V be an arbitrary, but fixed, vertex. We begin by proving that 
(A, M, v)->(B, z, z) for every (A, m, v) e KiG). 

Suppose, then, that (A, m, v) e HiG). We can generate a new state from (A, m, v) via the map / ; 
HiG) — > HiG) with /(A, M, v) defined by the following prescription: 
if du{A) - 1 then 

Choose an edge mm' € B with mm' ^ A 
return (A U mm', m', v) 
else if c/i (A) - 1 then 

Choose an edge vv' e B with vv' i A 
return (A U vv', u, v') 
else if <i„(A) = 2 and A B then 
Choose ww' e B with ww' i A 
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return (A U ww', w', w) 
else if c/„(A) = 2 and A = B then 

return (B,z,z) 
else if liu(A) = 3 then 

Choose the edge uu' i B 

return (A \ mm', m', v) 
end if 

The key observation to make is that (A, m, v)~^f{A, u, v) for every (A, m, v) e "RCG). Indeed, 
if du{A) 2 we simply have Pg „[(A,m, v) — > /(A,m, v)] > 0. Suppose instead that du(A) - 2, 
which impHes u - v. Lemma 1275] shows that (A, m, m)^^(A, w, w) for all m, w e y, so if A = B 
then we clearly have (A, m, v) = (B, u, u)<^{B,z,z) = /(A, m, v). On the other hand, if A B, 
then there must be at least one edge ww' which is in B but not in A. Again, Lemma |231 shows 
that (A, M, m)«~^(A, w, w), and in addition we have 

PcAi^' w) ^ (A U ww', w', w)] > 0, 

so that (A, w, w)-^(A U ww', w', w), and consequently (A, m, m)-w(A U ww', w', w). Therefore, we 
indeed have {A, u,v)^f{A, u, v), and in fact (A, m, v)~^/"(A, m, v) for any n e N, where /" - 
/ o / o ■ ■ ■ o / denotes n-fold composition of / with itself. 

Now, whenever A B, the state /(A,m, v) has either one more occupied B-edge, or one 
less occupied non-B-edge, compared to (A, m, v). Therefore, since there are only a finite num- 
ber of edges in G, if we start in any (A, m, v) e HiG) and apply / repeatedly, then we must 
eventually have /"(A, m, v) - (B,z,z), with n necessarily ^n/fe. It then immediately follows that 
(A, M, v)^(B,z, z). Finally, the reversibility of P^n now implies that in fact (A, m, v)<-^(B,z,z), 
and since this holds for all (A, m, v) e 'R(G) this implies f c,n is ergodic. □ 

Lemma 2.5. Let G be a finite, connected, bipartite cubic graph. For every A € 'F(G) and every 
pair x,y e V we have (A, x, x)^^(A,y,y) under Pq „. 

Proof. We begin by showing that (A, x, x)-^(A, x' , x') for all x e V and x' ~ x, where x' ~ x 
denotes that x' and x are adjacent (i.e. they are neighbors). Firstly, note that if xx' i A, then we 
simply have Pc,n[(A, x, x) — > (A U xx',x',x)] > and Pc,n[(A U xx',x',x) — > (A,x',x')] > 0, 
which immediately implies (A, x, x)~^(A, x', x') in this case. 

Therefore, let us consider x" ~ x with xx" e A. Lemma l276l guarantees that there exists an 
alternating path, f - ziZi - ■ ■ Zik, such that zi = x, zit - x", z; z,+i i A for / odd, and z; Zi+i e A 
for / even. The key observation is that it is always possible, via transitions of Pc,n, to move 
one defect along such a path while leaving the other defect fixed, which can be seen as follows. 
Since (A,zi,zi) e S, we can always make the transition (A,zi,zi) — > (A U z\Z2,Z2,Zi) when 
ZiZi i A. Furthermore, since |A| - \V\ and Z3 + x, it follows that (A U z\Zi,Zi,Z\) G 2) U and 
(A U z\Zi \ ZiZt„Zt„Z\) e T ■ In general, if the position of the first defect is zi with i even, then 
the corresponding state will be in £) U 0, and the transition that moves the first defect from z, 
to Zi+\ by deleting the edge will occur with strictly positive probability. Conversely, if the 
position of the first defect is zi with / > 1 odd, then the corresponding state will be in T, and 
since z,z,+i ^ A, the transition that moves the first defect from zi to z,+i by adding the edge ziZm 
will also occur with strictly positive probability. Consequently, 

(A,zi,zi)^(AA£(nz2*,zi) e£)U0. 
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Now we can move the second defect along !P, leaving the first defect fixed at zik- This has 
the eff'ect of flipping each edge back to its original state, so we arrive at (A, x", x"). Indeed, an 
analogous argument to that above shows that 

{Am'P\z2k,zx)-^{A^E{'P)AE{'P\z2k,Z2k) = (A,x",x"). 

Therefore, (A, x, x)~^{A, x' , x') for all x' ~ x. However, since this was proved for general x, it 
follows that precisely the same argument could be applied to show that (A, x', x')-^(A, x, x), and 
so we have (A, x, x)*^(A, x', x'). Finally, since G is connected, transitivity immediately implies 
that in fact (A, x, x)«-^(A, y, y) for every y e V. □ 

Remark 2.7. We remark that the dynamics of the worm algorithm along alternating paths dis- 
cussed in Lemma 12.51 is very similar to the use of alternating paths by graph theorists in the 
construction of maximal matchings; see e.g. IstIi . 

Lemma 2.6. Let G be a finite, connected, bipartite cubic graph. For every A € 'F(G), every 
x € V, and every x' ~ x, there exists a path Z\Z2 ■ ■ - Zik in G such that zi — x, Z2k — x' , Zi Zi+\ i A 
for i odd, and Zi Zi+i e A /or / even. 

Proof. We begin by noting that any path P between x and x' ~ x must have odd length, because 
P + xx' is a cycle and G is bipartite. 

Now, let A e 'FiG) and x e V, and suppose xx' i A, xx",xx"' e A. The path xx' is then 
trivially a path of the above form, with k - 1, so let us focus on constructing a path from x to x". 

Let ziZ2 ■ ■ ■ Z2k Z\ be a cycle in (V, A). Since it contains an even number of edges, we can 
colour half of them blue, and half of them red, in such a way that each vertex Zi is incident 
to precisely one blue edge and one red edge. For example, we can colour each edge 
blue if / is even and red if i is odd. In this way the edges alternate red, blue, . . . , red, blue as we 
traverse the cycle. Since the cycles in (V, A) are vertex disjoint, such colourings can be performed 
independently for each cycle. 

After performing such a colouring, each vertex in (V, A) is incident to precisely one red, one 
blue, and one vacant edge. If we now colour each vacant edge green, then we obtain a proper 
3 -edge-colouring of G. Suppose we now interpret the red edges as vacant, and the blue and green 
edges as occupied. Each vertex will again have degree 2 (one blue edge plus one green edge), so 
this procedure generates a new bond configuration Ared 6 'F(G). Furthermore, since each vertex 
is incident to precisely one green edge and one blue edge, each cycle ziZ2 ■ ■ ■ Z2k Z\ in Ared must 
be such that the edges alternate green, blue, . . . , green, blue as we traverse the cycle. 

Now, fix A e TiG) and x e V, and suppose xx' ^ A, xx" , xx'" e A. Suppose we perform a 
colouring of A, as described above, in which the edge xx" is blue. It then follows that there is 
a cycle z\Z2 ■ ■ ■Z2kZ\ in (y,Ai-ed) in which z\ - x, Z2 - x' and Z2k - x" ■ But this defines a path 
Z\Z2 ■ ■ ■ Z2k in G in which the edges are green (vacant in (V, A)) when / is odd, and blue 
(occupied in (V, A)) when / is even. We have therefore showed that there exists an alternating path 
of the required form between x and x" . A similar argument can obviously be used to construct 
the required path between x and x"'\ in fact the same colouring can be used as for the xx" path, 
provided we interpret the red and green edges as occupied, and the blue edges as vacant. □ 
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2.6. Connectivity-checking and the colouring method 

An important practical matter when implementing the algorithms we have presented so far, 
is the need, when n + 1, to perform a non-local query to determine if the cyclomatic number 
changes when an update is performed. Consider a spanning subgraph {Y,A) c G of a graph 
G — (y, E). Since the number of components, fe(A), is related to the cyclomatic number, c(A), by 
k(A) = |y| - |A| + c(A), the task of determining whether an edge-update changes the cyclomatic 
number is equivalent to determining whether it changes the number of connected components. 
The latter question can be answered by known dynamic connectivity-checking algorithms 115 If . 
which take polylogarithmic amortized time. A much simpler approach, which runs in polynomial 
time, but with a (known) small exponent is simultaneous breadth-first search lf52ll . We used the 
latter approach in the simulations presented in Section[3] In Sections [2.6.1l and l2.6.2l we discuss 
a different approach, the colouring method, which avoids altogether the need for such global 
queries, at least when n > 1 . Before discussing the colouring method, however, we make some 
remarks regarding the practical implementation of connectivity queries. 

To illustrate, we will consider Algorithm |2] with n > 1. Algorithmic states that connectivity 
queries are necessary only when (A, m, y) e U 2). In practice, however, even in this case one 
does not usually need to perform such queries. Suppose we assign a fixed (but arbitrary) ordered 
labeling to V, so V = {vi, V2, . . .), and let re [0, 1] be a uniformly-distributed random number. 
For notational convenience we set p = 1 /(n + 2) and e = 1 /3 - /? > 0. We can implement 
the du(A) - dy{A) = 3 block in Algorithmic as follows. Denote the neighbors of u by uj with 
j - 1, 2, 3, such that ui has the smallest label, and M3 the largest. If r e [{j - 1) p, j p] then we 
can simply choose uj without knowing the topology of the defect cluster. If instead r > 3 p, then 
we need to determine whether or not any of the edges uuj are bridges; there can be at most one. 
If MMj is a bridge then we choose uj, otherwise if none of the uuj are bridges then we choose uu j 
iff r G [3 p + (j - l)e,3 p + j e]. An analogous trick can be employed when n < 1. 

The computational burden imposed by these connectivity queries is of greater concern when 
jc < 00 than when x = 00. At x = +00, connectivity-checking is only required when du(A) - 
dy(A) = 3. By contrast, when jc < 00 it is in principle always necessary unless both defects are 
isolated, except where avoided by a trick of the type described above. For this reason, one might 
expect that the colouring method (whose raison d'etre is to avoid connectivity queries) would 
be more advantageous when x < oa. Our simulations suggest that this is indeed the case. In 
fact, while the colouring method was significantly more efficient on the critical branch, it was 
significantly less efficient than the connectivity-checking version when x - 00. 



2.6.1. The colouring method 

Consider a finite graph G = (V, E). The colouring method is a general method- 

ology for simulating models of the form 

0c,w(A)cx Y[ mC), AC£, (19) 

CeK(A) 

where K{A) denotes the set of all connected components of (V,A), and W is an arbitrary map 
that associates a nonnegative weight to every connected subgraph of G. Many lattice models in 
statistical mechanics can be expressed in the form ( fT9l ). If we set W{C) - q v'^*'-*! for all C, then 
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we recover the standard random-cluster model of Fortuin-Kasteleyn lfl9l] . If, instead, we set 



1, if C is an isolated vertex, 

nxl^^'^'l, if C is a cycle, (20) 
0, otherwise. 



and G has maximum degree A(G) < 3, then (pc,w = <pG,n,x and we recover the loop model ([T]). 

The key step in applying the colouring method is to choose an appropriate nonnegative weight 
function W < W for which we have a transition matrix ^ to simulate (f>Q^- In practice, 
appropriate means that ^ is more efficient/convenient to implement than any algorithm we 
have at hand for (pc.w- We shall return to this point in Section 12.6.21 Given W and PQ^,, the 
colouring method simulates (pc.w by an algorithm which, at each step, updates the bonds on a 
suitably-chosen random subgraph H c G using f ^ ^, while leaving all other bonds fixed. 

Algorithm 3 (colouring method). 
loop 

Current state is A Q E 

Independently colour each C e K{A) red with probability W{C)IW{C) and blue otherwise 
Identify the active subgraph Gred — G[yred] 
Choose a new A[^^ via P^^^^^^iA^^i -> A[^^] 
New state is A' — A'^^ U Abiue 
end loop 



By independently colouring each cluster in K(A) we obtain a random 2-colouring of the vertices 
cr 6 {red, blue)^ for which each edge in A has both its endpoints coloured the same colour 
We consider the subgraph induced by the red vertices Gied(o') - G[V'red(o')] as active and that 
induced by the blue vertices as frozen. The set A,ed = A n ^(Gied) is the subset of all edges in A 
for which both endpoints lie in Vied, and similarly Abiue -An E{Ghiue)- 

In Section l2.6.2l we specialize Algorithm[3]to construct a worm algorithm for the loop model 
</>c,n,x, which we present in Algorithm ID We conclude the current section by providing a more 
precise statement of Algorithm [3] in terms of transition matrices. To this end, let us set W - 
(Wredi Vl^biue) - {W, W - W) and iutroducc the following joint measure of colours and bonds 

/zc,w(A,£r)(xA(A,cr) ]~[ W^(c)(C), (21) 

CeK(A) 

where cr(C) denotes the colour of the vertices in cluster C e K(A), and A(A, cr) is the indicator 
for the event {(A,£r) : cr, = ctj for all ij e A). The transition matrix of the colouring method 
provides a Monte Carlo algorithm to simulate the joint measure ( 1211 1. Since the marginal measure 
2o- A'c,w(', f) on {A c E} is simply (pc.w, it then follows immediately that the colouring method 
in fact provides a Monte Carlo method for (pc.w- Proposition l2.7l provides a precise definition as 
well as a justification of the colouring-method transition matrix. Recall that the support of fic.w 
is, by definition, supp(yUG w) {iA,(r) : fic.'wiA,(r) > 0). In a slight abuse of notation, we also 
write cr e supp(yiG,vv) whenever there exists (A, cr) e supp(//G,w)- 

Proposition 2.7 (colouring Method). Consider W < W with supp{W) = supp(W), and for 
each cr e supp(iJ.c,w) let Pc j(o-) w transition matrix with state space [A^^A £ £(Gred(£''))) 
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and stationary distribution (p^ ^^^j ^, and suppose ^ is ergodic on supp{^^ ^). If we define 
transition matrices -Pcoiour and Phond on suppQjc.v/) by 



Pco^.uA(A,o■)^iA',o■')]=6A,A'HA,o■') U \ 

CeK(A) 

Pbond[(A, Cr) ^ (A', O-')] = Sa-,0-' A(A', 0-) 5am„„,A^,„ ^'c,,,«^),w>[^i-ed ^ A^J, 

f/zen P := Pcoiour ^'bond is ergodic and has stationary distribution i-ic,w- 

Proof. We simply sketch the proof; see for further details. Ergodicity follows by noting 

that there is a positive probability of consecutively colouring the whole graph red an arbitrary 
number of times, and then relying on the ergodicity of P(J^,■ Stationarity of yUc.w follows by 
observing that /ic,w is in fact stationary with respect to both Pcoioui- and Pbond separately, and these 
latter two facts can be easily verified by noting that Pcoiour[(A, c) — > (A', cr')] - 5a,a' /^cwCc'IA), 
and 

PC,W(A I 0-) = A(A, cr) <^c„_,(o-),w(^red) 0Gy^,^.(o-)_t(Ablue)- 

□ 

2.6.2. Applying the colouring method to worm algorithms for the loop model 

Since we do not need to perform connectivity checks when n - 1, we consider worm updates 



for the n - \ model to be convenient, and we know from experience 1134 12511 that they are also 
efficient. When n > 1 it is therefore natural to choose 

TV(C) = P*''' ifCisEulerian, ^^^^ 
10, otherwise. 

The resulting algorithm proceeds as follows: 

Algorithm 4 (coloured worm algorithm). 
loop 

Current state is A 

Colour each isolated vertex red 

Independently colour each loop red with probability l/n 
Identify G^ed 

Choose, uniformly at random, v e V{G^ed) 

Use n — I worm updates on Gred to make a transition (Ared, v, v) — > (A^^^, v', v') 
New state is A' — A'^^^ U Abiue 
end loop 

Finally, let us consider the « = 1 worm updates in a little more detail. Let H c G and consider 
the following transition matrix on C(//) 

PhAA ^ A'] -J; X ^'^."=i^-[^^' ^' ^) ^ ^')]' (23) 

v.y'eV 

where Pc.n. v is the restriction of (|9]l to the Eulerian subspace ( fTOt . Ph,x[A — > A'] is the probability 
that, starting in A e C(II), we pick, uniformly at random, a location for the defects, v, then 
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perform worm updates from (A, v, v) until we arrive at a new A' e C(H), regardless of the new 
location of the defects. It is clear that the row sums of Ph,x are correctly normalized, so that 
it defines a stochastic matrix, and it is also clear that Ph,x is in detailed balance with (pH,n=\,x- 
In addition, since PG,n,x is ergodic on any G when x < oo, it follows immediately that Ph,x is 
ergodic for all x < oo. The transition matrices i{a-)w required in Proposition 12.71 are then 
chosen to be PG„d(o-),.i- Analogous transition matrices can obviously be constructed from ( fTTT i 
and ( fT6b . For the fully-packed case, Ph,x=+oo is not necessarily ergodic on all possible subgraphs, 
however Proposition 12.71 only requires that it be ergodic for H - G, which is guaranteed by 
Proposition l2.4l 

Remark 2.8. Let us return to the case of the genuinely n-dependent connectivity-checking ver- 
sions of the worm dynamics, as discussed in Sections [2. II 12.31 and l2.4l We note that the right- 
hand side of ( |23] ). with n left arbitrary rather than fixed to 1, defines a perfectly valid alternative 
worm algorithm, in which an additional step is added: whenever the defects collide, uniformly at 
random choose a new vertex to move them both to. In particular, for the FPL model, if one were 
to add such a move then ergodicity could be proved without recourse to Lemmas 12.51 and 12.61 
However, we find empirically that there is no practical advantage to adding these moves, and we 
did not use them in our simulations (except, of course, when using the colouring method). 

3. Numerical results 

We simulated the loop model ([1]) on an L x L honeycomb lattice with periodic boundary 
conditions, using the algorithms described in Section |2l In particular, we used the genuinely 
n-dependent algorithms described in Sections 12.31 and 12.41 as well as the colouring algorithms 
described in Section 12.6.21 We shall refer to these two distinct versions as the connectivity- 
checking and colouring versions, respectively. 

We considered both the cases n < 2 and n > 2. The questions studied differed substantially 
in these two cases, since the exact phase diagram is known when n < 2, while no exact results 
are known at all for n > 2. For each choice of n, and each possible branch (when n < 2) we 
simulated at least seven (and up to eleven) different choices L, in the range 12 < L < L,nax- The 
values of L^ax used depended on the choices of x and n. They are summarized in Table[T] 



Table 1: Summaiy of dift'erent values of L„,nx used in each simulation. 



Critical 


Densely-packed 


Fully-packed 


;i > 2 


n 


0.5 1.0 1.5 2.0 
240 360 240 240 


0.5 1.0 1.5 
120 240 240 


0.1 1.0 1.25 1.5 1.75 2.0 
120 240 240 240 240 240 


3 10 
120 120 



For n > 2, we simulated at n = 3 and 10, with the aim of verifying the existence of a phase 
transition, the nature of its universality class, and obtaining accurate estimates of the critical 
points. For « < 2, we simulated on both branches of (HJi, as well as at x = +co. We focused 
on identifying the scaling exponents of five distinct classes of observables, characterizing loop 
lengths, face sizes, the magnetization of dual Ising spin configurations and its staggered analogue, 
and the return time to the Eulerian subspace. Although the latter observable is, by construction, 
defined on the full worm space S(G), rather than on the loop state space C(G), we shall see that 
it appears to have a deep connection with the loop model itself. 
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3.1. Ohservahles measured 



We measured the following observables in our simulations. All observables were measured 
only when the defects coincided, with the exception of the return time, 7", which is defined on 
the full worm chain. 

• The number of loops Ni{A) = c{A). 

• The number of bonds A/i, (A) = \A\. 

Note that on the fully-packed branch we trivially have Nb = \y\ = 2L^, so we only mea- 
sured this quantity on the critical and densely-packed branches. 

• The length of the largest loop Xi 

• The mean- square loop length 

I 

where the sum is over all loops /. 

• The size of the largest face Qi 

• The mean- square face size 

§2 := L-^ 2 \f? (25) 

/ 

where the sum is over all faces /. Every loop configuration A on the honeycomb lattice 
can be decomposed into a number of faces, each consisting of a collection of elementary 
hexagons, such that every pair of neighboring elementary hexagons which share an unoc- 
cupied edge in A belong to the same face. The size |/| of face / is then simply the number 
of elementary hexagons which it contains. 

• The dual Ising magnetization At. 

This is measured by assigning an Ising configuration to the dual triangular lattice in such 
a way that the loops on the honeycomb lattice form the domain boundaries of the Ising 
spin configuration. Such Ising configurations can be defined in an unambiguous manner 
whenever the loop configuration winds the torus an even number of times, and we therefore 
only measured this observable when the loop configuration was in this subspace of the 
cycle space. In such cases the spin configuration is unique (up to a global spin flip cr 1-^ 

-CT). 

• The sublattice dual Ising magnetization At,, for j = 1, 2, 3. 

Since the triangular lattice is tripartite, we can independently consider the Ising magneti- 
zation on each of its three sublattices. 

• The return time T to the Eulerian subspace C(G) x V. 

From these observables we estimated the following quantities: 

• The loop-number density n; := L~^{Ni) 
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• The loop-number fluctuation C/ :- L ^var{Ni) 



• The bond-number density on the critical and densely-packed branches rti, := L^^{Nh) 

• The bond-number fluctuation on the critical and densely-packed branches C,!, := L^^\ai(Nb) 

• The expectations (Xi) and 

• The expectations (^1) and {^2) 



• The dual Ising susceptibility ;t'ising = L ^ ^^ising^ 



The second and fourth moments, and (Al^tag)' of the the staggered dual Ising 

magnetization, Alstag, defined by 

A^Lg := (Ml - Mif + {Ml - M^f + (M3 - Mxf 
The staggered susceptibihty Xstn% - l^^^ (^Lg) 



'stag 

• The dimensionless ratio 

/ A,, 

stag 



<A^lg) 



• The mean return time to the Eulerian subspace (7"). We also estimated the higher-order 
moments (7~*), for A: = 2, 3, 4, and the distribution P(7" = 0- 

For each quantity Y = n,, n^,, Q, Ch, <Xi), {Li), {Q\), {Qi), Arising, ;tstag, Qs and (7"*) we 
performed a least-squares fit of our Monte Carlo data to the finite-size scaling (FSS) ansatz 

y08,L) = CO + C1O6 -Pc) + ■■■ + L'^ [ flo + ai(j3-j3,)U' + aiifi - p.fl^y + ... 

(lb) 

+ biL-'^' +b2L-'^^ + ...]. 

Here y, is the leading thermal exponent. The a, are coefficients of the FSS variable (J3 - l3c)L!", 
the b; are the coefficients of the corrections-to-scaling terms, and the c, are the coefficients of 
analytic terms. There are also cross-terms involving products of terms arising from each of these 
three sources. The choice of which terms to include in the fit for a given choice of observable 
varied from case to case, and involved a certain amount of trial and error. The exponent xy in 
( l26l l is a generic label for whatever the dominant exponent happens to be for the quantity Y. In 
particular, if Y happens to be dimensionless, such as Qs, then we have xy -0 identically, and in 
this case all the c, are identically zero. 

As a precaution against corrections to scaling, we imposed a lower cutoff L > L^in on the 
data points admitted to the fit, and we studied systematically the effects on the fit of varying the 
value of LmiiT We used the Levenberg-Marquardt algorithm to perform the fits. 

3.2. Fits for n < 2 

The results of the fits for n < 2 are presented in Tables |2l [3] and |4] We discuss these results 
observable by observable in the following sections. 
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Table 2: Critical exponents for the critical branch. The exponent denotes the scaling dimension of {T) for the 

connectivity-checking version of the worm dynamics. 



n 


-^loop 


^hull 




Xpj, 






-^energy 




Xpj2 


-^worm 


Xh 


0.5 


0.648(2) 


0.6478 


0.0258(4) 


0.0261 


0.0567(2) 


0.0567 


0.816(3) 


0.79(2) 


0.8178 


0.1157(3) 


0.1154 


1.0 


0.623(2) 


0.6250 


0.0518(4) 


0.0521 


0.1251(2) 


0.1250 


0.998(3) 


0.97(2) 


1.0000 


0.1250(2) 


0.1250 


1.5 


0.592(2) 


0.5935 


0.0800(4) 


0.0801 


0.2194(2) 


0.2195 


1.254(5) 


1.28(6) 


1.2519 


0.1322(2) 


0.1322 


2.0 


0.5000(2) 


0.5000 


0.1249(3) 


0.1250 


0.4998(3) 


0.5000 


1.996(7) 




2.0000 


0.1251(5) 


0.1250 



Table 3: Critical exponents for the densely-packed loop branch. The exponent X„orm denotes the scaling dimension of 
(7~) for the connectivity-checking version of the worm dynamics. No exponents are reported forxsag since it was found 
to be constant. 



n 






-^face 


Xpj, 


-Rising 


Xpj 


^worm 


Xi, 


0.5 
1.0 
1.5 


0.139(2) 

0.2500(1) 

0.3508(3) 


0.1386 
0.2500 
0.3506 


0.0640(4) 
0.1044(4) 
0.1282(3) 


0.0637 
0.1042 
0.1280 


1.58(1) 
0.947(3) 


1.5843 
0.9482 


-0.0788(4) 
0.0000(1) 
0.0620(3) 


-0.0791 
0.0000 
0.0619 



3.2.1. Scaling of energy-like quantities 

Standard finite-size scaling arguments predict that on the critical branch 

nb ~ a + bV'''"^''^, (27) 
Ch ~ a + bL^y^"'''-^, (28) 

where jenergy = 2 - Xgnergy IS the fractal dimension characterizing the energy-like quantities, 
and ^energy is the corresponding scaling dimension. It was observed in 12411 that for the critical 



loop model with 1 < n < 2, we have ^energy = Xp,2 where Xp,2 is the second thermal scaling 
dimension of the <7-state Potts model with q - n^, whose expression in terms of the Coulomb gas 
coupling g is known Ii26il to be 

Xp,2 = -2 + 16/g. (29) 

For n - 1, one has Xpa - 1 and 2yenergy -2 = 0, and Eq. (l28l l is replaced hy Cb ~ a + b In L. 
Table|2]shows the numerical results, which confirms Eq. ( l29l l. 

The behaviour of the loop-number density n; and its fluctuation C/ is also described by 
Eqs. (|27]| and (|28]|. This is perhaps unsurprising in light of the Euler relation c(A) - \A\-\V\+k(A). 

On any graph of maximum degree 3, the observable Nh/Ni is simply the arithmetic mean 
of the loop lengths in a given configuration. This quantity was studied in detail in [53J, and in 
particular it was found that its expectation {NblNf) tends to a constant in the thermodynamic 
limit. 



Table 4: Critical exponents on the fully-packed branch. The exponent denotes the scaling dimension of {T) for 

the connectivity-checking version of the worm dynamics. 



n 


-^loop 






Xpj, 




Xpj 


-'''.stag 




-^worm 


Xi, 


0.1 


0.031(3) 


0.0309 


0.0154(4) 


0.0152 


1.92(6) 


1.9074 


0.322(5) 


0.3231 


0.032(4) 


0.0309 


1.0 


0.2500(2) 


0.2500 


0.1042(3) 


0.1042 




1.2500 


0.2499(1) 


0.2500 


0.2498(4) 


0.2500 


1.25 


0.3005(2) 


0.3006 


0.1181(2) 


0.1180 


1.10(3) 


1.0982 


0.2333(3) 


0.2331 


0.3005(3) 


0.3006 


1.5 


0.3505(4) 


0.3506 


0.1276(4) 


0.1280 


0.93(2) 


0.9482 


0.2165(3) 


0.2165 


0.3507(3) 


0.3506 


1.75 


0.4039(3) 


0.4042 


0.1336(3) 


0.1335 


0.79(1) 


0.7875 


0.1986(4) 


0.1986 


0.4043(5) 


0.4042 


2.0 


0.482(3) 


0.5000 


0.1345(3) 


0.1250 


0.57(2) 


0.5000 


0.175(3) 


0.1667 


0.480(3) 


0.5000 
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On the critical branch, the fitting results for the constants, riho, Cbo, n/o, Qo, and B - n nho/njo, 
are shown in Table|5l Our estimate B(n = 2) = 38.832(2) agrees well with 38.834(2) in Ref. E^ . 

On the densely-packed branch, our data indicate that the scaling of «/ and C; fits well the 
formula rii ~ a + bL^^ and Ci ~ a + bL^^, while for rit, and C/,, there's no detectable finite-size 
dependence. The asymptotic behaviour of B for n — > agrees with the prediction B{n — » 0) = 
35.70(2) m. 

On the fully-packed branch, we have Nt, = \V\ and so the expectation of the arithmetic mean 
of the loop lengths is simply 2/n;. Our data for «/ fit the formula rii ~ a + bL^^, however, for 
n - 2 the dominant correction exponent appears to be -2.20(4), most probably arising from 
logarithmic corrections. No detectable finite-size dependence is found for C/. For n - 2, the 
estimates for n;o, C/o and B confirm the prediction «/ = 1/9, Ci - 1/9 + 1/135, and B = 36. We 
also determined the asymptotic value B(n — > 0) = 30.04(2), in good agreement with the exact 
value 30.0344 . . .. 



Table 5: Numerically determined results fij,o, Ctio, n;o, C;o, and B = tmijo/nio. 



Branch 


n 


ni,o 


Cbo 


n/o 


Cio 


B 




0.5 


0.26646(6) 


-2.37(4) 


0.018644(5) 


0.011(4) 


7.146(4) 


Critical 


LO 


0.49998(2) 


0.41(4) 


0.033664(4) 


-0.0072(3) 


14.852(3) 




L5 


0.72953(1) 


4.8(1) 


0.046243(3) 


0.078(2) 


23.664(3) 




2.0 


1.114837(3) 


1.789(2) 


0.057418(2) 


0.0553(2) 


38.832(2) 




L5 


1.395506(2) 


0.9760(3) 


0.0475624(4) 


0.03978(2) 


44.0108(6) 


DPL 


LO 


1.499999(2) 


0.7500(1) 


0.0352504(4) 


0.02999(2) 


42.5527(6) 




0.5 


1.579679(2) 


0.5960(2) 


0.0199245(6) 


0.01798(2) 


39.642(2) 




0.1 


2 





0.006527(4) 


0.00640(2) 


30.642(9) 




LO 


2 





0.057668(2) 


0.05242(3) 


34.6813(4) 


FPL 


1.25 


2 





0.070680(1) 


0.06484(4) 


35.3707(2) 




1.5 


2 





0.083677(2) 


0.07851(5) 


35.8521(4) 




1.75 


2 





0.096989(2) 


0.0952(2) 


36.0866(4) 




2.0 


2 





0.111111(2) 


0.1185(1) 


36.0000(3) 



3.2.2. Scaling o/(Xi) and {£,2) 

Standard finite-size scaling arguments predict that 

(£,)^LV,o„p (30) 
(£2) ~ L^-""-"^ (31) 
where yioop - 2 - X\oo^p is the fractal dimension characterizing loop length, and Xioop is the 



corresponding scaling dimension. It was argued in Il54ll that for the critical loop model we have 



^loop - ^huih where Xhuii is the hull scaling dimension 

^huii = 1 - (32) 

g 

and g is related to x and n as in (|5]l and (|6]l. It was argued in [46] that X\oo^ - Xhuii also holds on 
the fully-packed branch. 

On the critical branch, Xioop = ^huU was verified for 1 < « < 2 in the Monte Carlo study 
presented in f24]. Table |2] confirms this result, and shows that it extends to n < 1. Since the 
expression (l32l i for Xhuii(g) as a function of g should be universal, we expect that if we insert the 
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DPL expression for g into ( l32l l then we would again have X\oop - ^huii- Table [3] shows that this 
is indeed the case, and Table |4] shows that the result also holds on the fully-packed branch. In 
Fig.[3]we plot our numerical estimates of Xioop together with the exact result for Xhuu for all three 
branches. The agreement is clearly excellent. The small deviation of the FPL estimate aXn - 2 
is presumably due to the presence of logarithmic corrections to scaling. 

3.2.3. Scaling of{@\) and {Qi) 

0.14 
0.12 
0.1 

^ 0.08 
0.06 
0.04 
0.02 





0.5 I 1.5 2 

n 

Figure 4: Numerically determined scaling dimension Xf^cc plotted with the exact expression )35t for Xpj,, as a function 
of n. 




Predicted 
Numerical 
Numerical on FPL 



Analogously to the previous case for loop length, we expect that 

{@x)~V^-' (33) 
{Qi) ~ (34) 
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where yf^^e = 2 - Xf^^e is the fractal dimension characterizing face size, and Xface is the corre- 
sponding scahng dimension. It was argued in Ii55.1 that for both the critical and densely-packed 
loop models we have Xf^ce - Xpj,, where Xp^i, is the magnetic scaling dimension of the g'-state 
Potts model with q - n^, whose expression in terms of the Coulomb gas coupling g is known f2^ 
to be 

{6-g){g-2) 



X, 



P,h 



(35) 



and g is related to x and n as in (|5]l and (|6]l. 

On the critical branch, Xface - Xpj, was verified for 1 < n < 2 in the Monte Carlo study 
presented in (l^. Table |2] confirms this result, and shows that it extends to « < 1, and Tabled 
verifies it in the densely-packed phase. It is not clear a priori that this relationship should also 
hold in the fully-packed branch, but Table |4] shows that it does. In Fig.|4]we plot our numerical 
estimates of Xf^ce together with the exact result for Xp^i, for all three branches. The agreement is 
clearly excellent. The deviation of the FPL estimate at n = 2 is presumably due to the presence 
of logarithmic corrections-to-scaling. 



3.2.4. Scaling of ximng 




Figure 5: Numerically determined scaling dimension Xisin„ plotted with the exact expression i37\ for Xpj, as a function 
of II. 



It is natural to expect that 

;risi„g ~ L^'^""'""' (36) 

for some scaling dimension X^^ng- It is not a priori obvious what the form of Xjsing should be, 
however it was observed in 12411 for 1 < « < 2 that on the critical branch we have X^ing = Xp^,, 
where Xp, is the thermal scaling dimension of the ^-state Potts model with q = n^, whose 
expression in terms of the Coulomb gas coupling g is known L26.1 to be 



Xp,^--l, 



(37) 



and where g € [4, 6] is related to x and n as in (|5]l. Table |2] confirms this result, and shows that 
it extends to n < 1 . One would expect that the relationship would continue into the densely- 
packed phase, with g e [2,4], and Table [3] shows that this is indeed the case. It is not entirely 
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obvious what the behaviour should be on the fully-packed branch, but Table |4] shows that it 
simply coincides with that of the densely-packed branch. We note that on the densely-packed 
and fully-packed branches, when x < V2 we have 2 - 2Xising < and so in these cases we 
estimated Xising from an ansatz of the form 

Xung^a + bL^-^''''-'^ + ... (38) 

in which the i}-^^i^m term is sub-dominant. 

In Fig. |5]we plot our numerical estimates of Xising together with the exact result for Xp, for 
all three branches. The agreement is clearly excellent. We note that there is no data point for 
n = 1 on the densely -packed branch, since (x, n) - (1,1) simply corresponds to site percolation 
on the dual triangular lattice (the +(-) spins are regarded as occupied (vacant) sites). In this case, 
the nontrivial dependence of;t'ising on the system size vanishes; i.e. the amplitude gq in ( |26] ) is 
identically zero. 

3.2.5. Scaling of X slug 

0.34 I , , , 1 




0.16 I ' ' ' 1 

0.5 1 1.5 2 

n 

Figure 6: Numerically determined scaling dimension Xstaa as a function of n in the fully-packed phase, plotted with the 
conjectured exact expression Xstag = 2/3g. 

On the critical branch, one might intuitively expect that x&vdg converges to a constant as L 
increases, since the symmetry of the sublattices should cause the staggered magnetization to 
cancel out for ferromagnetic models. Indeed, our numerical data show that on the critical branch 
Afstag is well described by the simple FSS ansatz 

X.t,.^a + bL-''^^'K (39) 

Perhaps surprisingly, our data also strongly suggest that we can in fact make the identification 
^stag = ^p,t2 where Xpa is the second thermal scaling dimension of the ^-state Potts model with 
q - n^, whose expression in terms of the Coulomb gas coupling g is ( |29l ). See Table|2l 

In the densely-packed phase, when n - 1 it is clear that stag is simply a constant which dis- 
plays no finite-size dependence, since the model corresponds to site percolation on the triangular 
lattice in this case. It is not obvious that this behaviour should persist when n ^ 1, however 
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our simulations strongly suggest that this is indeed the case: we find empirically that there is no 
finite-size dependence ofxstug for any < n < 2 on the densely -packed branch. 

In the fully-packed phase, the n = 1 case corresponds to the zero-temperature antiferromag- 
netic Ising model on the dual triangular lattice, which is known to be critical, and;i'stag is its order 
parameter. It follows that for n = 1 we have 

X.i.s ~ L^'^""'"' (40) 

with Xstag = 1/4. Our numerical results show convincingly that in fact ( l40l l holds for all n < 2. 
However, we were not able to identify the exponent Xst^g in terms of other known exponents. 
Instead, we fitted our empirical estimates of Xjtag to a simple Coulomb-gas ansatz 

Xis) = ^^^^ (41) 
dg 

with fl, b, c, d unknown integers. We found that Xstag is well described by the simple formula (see 
Fig.© 

^stag = (42) 

Based on the agreement presented in Table [3] we conjecture that this formula is in fact exact. 
3.2.6. Scaling of {T) 

0.5 
0.4 
0.3 

0.2 

0.1 

-0.1 
-0.2 

0.5 1 1.5 2 

n 

Figure 7: Numerically determined scaling dimension X„o,,n plotted with the exact prediction for X;, as a function of n. 




It is natural to expect that 

{T) ~ L2>v„,m-2 (43) 

where jworm = 2 - X„orm is a fractal dimension characterizing the time of return to the Eulerian 
subspace C(G) x V c S{G) in the connectivity-checking version of the worm dynamics. 

Although T is inherently an observable on the state space of the worm dynamics, S{G), 
rather than on the state space of the loop model, C(G), its scaling is linked to the loop model in 
a fundamental way. In fact, we find numerically that Xv/om - Xf, for all three branches and for 
all n < 2, where Xk is the magnetic scaling dimension of the 0(n) loop model. We note that, as 
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shown in on the critical branch of the n - I model {7") is equal to the susceptibility of 

the Ising model whose high-temperature graphs are sampled by the worm dynamics. Therefore, 
in this special case we have a sound theoretical argument that X„orm - Xi„ so the general result 
is not completely unexpected. 



On the critical and densely-packed branches, the exact expression 11431 14111 for X/, is given by 



2 3 

Xh^l —g, (44) 

g 32 



where g is related to x and n as in © and while on the fully-packed branch it is i44ll45ll46ll 



Xh^l- (45) 

8 

with g e [2, 4] related to n as in (|5]). Tables|2][3]and|4]show a comparison between Xwoim and Xh, 
and in Fig. [T] we plot our numerical estimates of X„orm together with the exact result for Xi, for 
all three branches. The agreement is clearly excellent. 

Based on these results for the honeycomb lattice, it appears that the identity Xworm - X/, is 
a generic property of the worm algorithms we present. Preliminary results for the hydrogen- 
peroxide lattice further support this conjecture; these results will be reported in detail elsewhere. 
For this reason, it seems in a certain sense that the worm algorithm provides the natural dy- 
namics for the 0(n) loop model. The relation Xworm - X/, also has some important practical 
consequences, since it provides a simple way of directly measuring X/,. As one application, it 
can be used to obtain an accurate estimate of the curve X/,(n) as a continuous function of n in 
three dimensions, where no exact results are currently known. 

The identification X„orin - Xi, also has implications for the efficiency of the worm dynamics. 
For n < 1, the exponent 2 - 2Xh governing (7~) is greater than 2 in the densely -packed phase, 
implying that it takes a time larger than order volume between visits to the Eulerian subspace. 
The efficiency of the worm dynamics will therefore suffer in this region. By contrast, on the 
critical branch one has 2 - IX/, e (1 .73, 1 .8) for all < n < 2, while on the fully-packed branch 
2 - 2Xh decreases monotonically with n e [0, 2] from 2 to 1 . 

Remark 3.1. Although we have not proved that Algorithm|2]remains ergodic at « = 0, we used it 
to perform simulations of the uniform Hamiltonian cycle model, and we found that X„orm = = 
Xi,(n = 0). This provides some modest evidence that the conjecture Xwoim - Xi, in fact extends 
to n = when jc = oo, and also that Algorithm |2] may indeed be ergodic. 

Remark 3.2. For both the critical branch and densely -packed branch, one has the choice of 
using either the simple worm dynamics, given by (|9]l, or the rejection-free worm dynamics ( fTTT l. 
Perhaps unsurprisingly, both versions have the same value of Xworm- On the fully -packed branch, 
only the rejection-free dynamics can be used, since the simple dynamics is non-ergodic. 

Remark 3.3. When n > 1, one can also consider Xwomi for the colouring version of the worm 
dynamics. On the critical branch, we simulated both versions and measured T in both cases. 
There appears to be no reason, a priori, for (7~) in the two cases to be related in any simple way, 
or even to share the same value of Xworm- However, empirically, we found that on the critical 
branch the exponent X„orm governing (T) agrees, within error bars, for the two versions. On 
the densely-packed branch, by contrast, we found Xwoim - 0.053(2) for the colouring version 
when n - 1.5, which appears to disagree with the value X„orm - 0.0620(3) for the connectivity- 
checking version reported in Table [3] Note that this implies the colouring version has the larger 
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Table 6: Estimated values of the scaling exponents, y and t, appearing in the ansatz j46t for the distribution P(7~ = t). 
Values for the critical, densely-packed, and fully-packed branches are shown. 



Critical 


Densely-packed 


Fully-packed 


n 


y yworm 


^ y Vworm 


^ y yworm 


0.1 
0.5 
1 

1.5 

2 


1.286(4) 2.479(4) 1.8843(3) 
1.248(3) 2.325(3) 1.8750(2) 
1.223(4) 2.234(4) 1.8678(2) 
1.174(4) 2.120(4) 1.8749(5) 


1.044(5) 2.257(7) 2.0788(4) 
1.077(4) 2.164(6) 2.0000(1) 
1.108(5) 2.103(8) 1.9380(3) 


1.166(6) 2.32(1) 1.968(4) 

1.364(3) 2.360(5) 1.7502(4) 
1.471(4) 2.458(7) 1.6493(3) 
1.602(5) 2.62(1) 1.520(3) 



mean return time (7~), since 2 - 2X„orm is larger in this case. Since we only simulated one value 
of n > 1 on the DP branch, this is the only point where a comparison between the two versions 
can be made. We did not perform any systematic numerical tests of the colouring algorithm in 
the fully-packed phase, since our preliminary results suggested it was significantly less efficient 
than the connectivity-checking version. 

3.2.7. Scaling of V{T = f) 

In order to better understand the behaviour of T, we computed its histogram and thereby 
estimated its distribution ¥{T = t). One expects that F(T - t) obeys a scaling law for large t in 
a critical system. More precisely, we expect that at finite volume we have 

P(r = f)~r7(f/L0, (46) 

for some choice of exponents t and y, and some scaling function /. Note that the ansatz (|46] | 
implies that (7~*^) ~ to leading order as L ^ oo. 

The distribution of the loop length was studied in |46] and found to be of the form (l46b . We 
have computed the distributions of X2 and Q2 in the present work, and also found excellent fits to 
( l46l l. In addition, distributions of the form ( |46] ) are known to govern the cluster-size of the random 
cluster model [52]. For all these cases, the exponent y appearing in ( |46] | coincides with the fractal 
dimension of the corresponding geometric objects; in particular we have y - yace = 2 - Xph and 
y - yioop = 2 - Xhuii for the face-size and loop-size cases, respectively. In such cases the two 
parameters y and t appearing in ( |46] ) can be related by the following argument. Consider a 
generic observable characterizing a geometric property (cluster-size, loop-size, face-size. . . ) 
with fractal dimension yy\. Since a given object has scale U^, the probability of picking one at 
random should be ~ /U' so that we expect {Jl'') ~ U^'^ L*^-*'" = /^(*+i).v.^-2 Combining this 
latter result with the fact that (^1*) ~ iji<+'^-T)y ^ ^jj^j assuming y - yj\, it follows that r - Ijy. 

Table |6] lists our fits for y and r, together with ywoim for ease of comparison, and Fig. [8] 
shows a finite-size scaling plot of P(7" > t) for n = 1 on the critical branch. Contrary to the 
observables discussed above, our numerical results show clearly that for the observable T, the 
exponent y appearing in ( |46] | is not equal to the corresponding fractal dimension y„orm- Although 
the combined exponent (2 - T)y -2 - 2Xi, governing the mean (T) is universal, it remains an 
open question whether or not t and y are themselves universal. We tried to fit the values of y in 
Table|6]to the ansatz (HTI) with d - 8, 12, and 16, but we did not obtain a meaningful expression 
within the estimated statistical errors. 

3.3. 0{n) loop model for n > 2 

Until quite recently, it was generally believed lf56|] that the (two-dimensional) 0{n) loop 
model, ([TJ, does not exhibit a phase transition when n > 2. Using numerical transfer matrix 

28 




10"' 10"^ 10'' 10° lO' 10^ 



Figure 8: Finite-size-scaling plot showing f¥{T > t) versus t/D' for « = 1 on the critical branch, with system sizes 
L = 12, 24, 48, 96, 192, 384, 768. 



Table 7: Estimates of X,, Xj, and for the 0(n) loop model with n = 3, 10. For comparison, note that for the critical 
3-state Potts we have X, = 4/5 and Xi, = 2/15. 



n 


X, X, Xi, 


3 

10 


6.822(7) 0.81(3) 0.132(3) 
1.5430(2) 0.83(4) 0.135(4) 



methods, however, strong evidence was found in |21|] that a line of critical points does in fact exist 
for all n > 2. Furthermore, the results presented in 12111 also suggest that this phase transition 
falls into the three-state Potts universality class. This can be understood by noting that the n — > oo 
loop model is equivalent to the hard-hexagon model |7]. 

We performed Monte Carlo simulations of the n = 3 and « = 10 loop models, using the 
worm algorithm presented in Section |2l We computed the staggered magnetization on the dual 
triangular lattice, which allowed us to then compute the staggered susceptibility and the dimen- 
sionless ratio Qs- By studying Qs we could accurately locate the critical point. See Fig. |9] and 
Table |7] By studying both stag and Qs we also estimated the exponents X, and X/,. See Table |7] 
Our numerical values for X, and X/, are entirely consistent with a transition in the 3-state Potts 
universality class, which has X, = 4/5 = 0.8000 and Xi, = 2/15 = 0.1333. 

For comparison, we note that in [21] it was estimated that Xc = 1.52(1), X, = 0.80(1) and 
Xh - 0.14(1) for n - 10. No data for n - 3 was reported in |2lll, however their quoted values 
for n = 4 were X, = 0.76(5) and X/, - 0.1(1). A number of different n values were reported 
in 112 ill , and the accuracy of their estimated exponents was found to increase with n. By contrast, 
we found the computational effort required for the n - 3 and n = 10 cases to be comparable. 



4. Discussion 

We have presented a Markov-chain Monte Carlo algorithm of worm type that correctly simu- 
lates the 0(n) loop model on any bipartite cubic graph, for any n e (0, oo) and x e (0, oo], and we 
have proved rigorously that the algorithm is ergodic and has the correct stationary distribution. 
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We have then applied this algorithm to the honeycomb-lattice loop model. Comparing our nu- 
merical results when n < 2 with Coulomb gas theory allowed us to identify the exact exponents 
of a number of fundamental geometric observables, as well as observables related to dual Ising 
spins. Furthermore, we have provided compelling numerical evidence that Xworm - X/, in all three 
branches for all n < 2. This suggests that Xwoim can be used as an efficient means to estimate X/, 
in models where no theoretical predictions exist, such as in three dimensions. For the case n > 2, 
we confirmed the existence of a phase transition in the 3-state Potts universality class, which has 
previously been observed using transfer matrices. Equipped with our worm algorithm, it is now 
natural to repeat these studies on bipartite cubic graphs other than the honeycomb-lattice, includ- 
ing to the Hydrogen-peroxide lattice, which provides a natural three-dimensional generalization, 
and the (4 ■ 8^) Archimedean lattice, which is the dual of the Union Jack lattice. These results 
will be reported elsewhere. 

Another natural application of the worm algorithms presented in this article is to the study of 
antiferromagnetic Potts models on triangulations of the torus. Since antiferromagnetic models 
do not display universality, it is of significant interest to study such models on a variety of dif- 
ferent lattices. As discussed in the introduction, it is well known that the honeycomb-lattice FPL 
model with n = 1 is equivalent to the zero-temperature triangular-lattice antiferromagnetic Ising 
model. While cluster algorithms iEtIEsIi for this model are thought to be non-ergodic, the worm 
algorithm presented in Section|2]provides a provably valid Monte Carlo method; see [25]. 

This observation can be generalized in two ways. Firstly, the dual of any bipartite cubic map 
on the torus is an Eulerian triangulation, and the worm algorithm from Section |2] can immedi- 
ately be applied to study Ising models on these triangulations. Secondly, as already mentioned 
in the Introduction, the honeycomb-lattice FPL model with n = 2 is equivalent to the zero- 
temperature triangular-lattice 4-state Potts antiferromagnet ll29ll. as well as the zero-temperature 
kagome-lattice 3-state Potts antiferromagnet. As also already mentioned, the Wang-Swendsen- 
Kotecky fsfl] (WSK) cluster algorithm, which is undoubtedly the current state-of-the-art for sim- 
ulating antiferromagnetic Potts models, has recently been proved iH , 32 ] to be non-ergodic for 
both of these models. By contrast, the worm algorithms described in Section|2]have been proved 
to be ergodic, and they can be easily applied to the corresponding loop models, in order to then 
study these Potts antiferromagnets. To do so, it is sufficient to simply augment the usual worm 
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steps with an additional step that randomly assigns an alternating colouring (red, blue, red,. . .) 
to the edges of each cycle. Each cycle can be coloured in precisely two ways. By interpreting the 
remaining edges as green, this defines anew transition matrix on 3-edge colourings of the honey- 
comb lattice, given by 2~'^^'*'f c,n=2[(^, m, v) (Ahuu' , u' , v)]. Since f g,« is in detailed balance 
with (pcni-A) oc n^*'*', at « = 2 this new transition matrix uniformly samples 3-edge colourings. 
But since the kagome lattice is the medial graph of the honeycomb lattice, it immediately fol- 
lows that this new transition matrix in fact uniformly samples 3-vertex colourings of the kagome 
lattice (i.e. simulates the zero-temperature kagome-lattice Potts antiferromagnet). The 4-state 
model can be treated in a similar way. Finally, we note that these mappings from n = 2 loop 
models to zero-temperature q - 3 and q - 4 state Potts models hold quite generally. In particu- 
lar, the 4-state antiferromagnetic Potts model on a variety of two-dimensional triangulations can 
be studied using such worm algorithms. AH of these possibiUties we leave for future work. 
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